load packages

library(tidyverse)
library(knitr)

define palettes

pal_control = c("#2E86B4", "#F2B827")
pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_condition = wesanderson::wes_palette("Zissou1", 6, "continuous")
pal_condition8 = wesanderson::wes_palette("Zissou1", 16, "continuous")
pal_condition8 = c(pal_condition8[1], pal_condition8[3],
                  pal_condition8[14], pal_condition8[5],
                  pal_condition8[15], pal_condition8[16],
                  pal_condition8[9], pal_condition8[12])

source data

Note that % fat variable was incorrect in 2012_FDES –> replaced it with the data from the EMA dataframe

source("load_data.R")

ind_diffs = ind_diffs %>%
  mutate(bmi_z = as.numeric(scale(bmi)),
         fat_z = as.numeric(scale(fat)))

standardize brain data

Windsorizing observations +/- 3 SD from the mean to 3 –> otherwise we lose these subs in the SCA analyses
Z-score within ROI and session for ROIs
Z-score within map, test (association or uniformity), mask (masked or unmasked) and session for PEVs

# roi distributions
betas %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, roi, process, meanPE) %>%
  unique() %>%
  filter(!is.na(roi)) %>%
  ggplot(aes(meanPE, fill = roi)) +
  geom_density(color = NA) +
  facet_wrap(~roi, ncol = 4) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

# dot product distribution
dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = test)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct, mask) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = mask, alpha = .5)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

# standardize and winsorize
betas_std = betas %>%
  group_by(roi, session) %>%
  mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
         outlier = ifelse(meanPE_std > 3, "yes",
                   ifelse(meanPE_std < -3, "yes", "no")),
         meanPE_std = ifelse(meanPE_std > 3, 3,
                      ifelse(meanPE_std < -3, -3, meanPE_std))) %>%
  ungroup()

dots_std = dots %>%
  group_by(map, test, mask, session) %>%
  mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
         outlier = ifelse(dotProduct_std > 3, "yes",
                   ifelse(dotProduct_std < -3, "yes", "no")),
         dotProduct_std = ifelse(dotProduct_std > 3, 3,
                          ifelse(dotProduct_std < -3, -3, dotProduct_std))) %>%
  ungroup()

summarize outliers

# summarize
betas_std %>%
  filter(session == "all") %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))

dots_std %>%
  filter(session == "all") %>%
  group_by(outlier, test, mask) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))

# remove outlier column
betas_std = betas_std %>%
  select(-outlier)

dots_std = dots_std %>%
  select(-outlier)

merge data

dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
  mutate(subjectID = as.character(subjectID)) %>%
  left_join(., ind_diffs, by = "subjectID") %>%
  ungroup() %>%
  mutate(subjectID = as.factor(subjectID),
         condition = as.factor(condition),
         control = as.factor(control),
         roi = as.factor(roi),
         process = as.factor(process),
         test = as.factor(test))

saveRDS(dataset, "full_dataset.RDS")
saveRDS(betas_std, "betas_dataset.RDS")
saveRDS(dots_std, "dots_dataset.RDS")

tidy data for plotting

Select cons generated across all sessions
Calculate mean across ROIs within each psychological process
Calculate balance score ROI = cognitive control - reward

dots_plot = dots_std %>%
  filter(mask == "unmasked") %>%
  filter(session == "all") %>%
  select(-c(map, con, mask, session, dotProduct)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender))

betas_plot = betas_std %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE)) %>%
  select(-c(con, xyz, roi, meanPE, sdPE, meanPE_std)) %>%
  unique() %>%
  spread(process, meanProcessPEstd) %>%
  mutate(balance = cognitive_control - reward) %>%
  gather(process, meanProcessPEstd, cognitive_control, reward, value, balance) %>%
  filter(session == "all") %>%
  select(-c(session)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  ungroup()

combined_plot = dots_plot %>%
  rename("value_std" = dotProduct_std) %>%
  mutate(test = ifelse(process == "craving_regulation", "neural signature", test)) %>%
  select(subjectID, process, test, condition, control, value_std, fat, bmi) %>%
  bind_rows(betas_plot %>%
              rename("value_std" = meanProcessPEstd) %>%
              mutate(test = "ROI") %>%
              select(-contains("_z")))

tidy data for modeling

ema_tidy = ema %>%
  select(subjectID, desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m, hunger.m) %>%
  unique() %>%
  mutate(enact_prop_z = as.numeric(scale(enact_prop)))

model_data = dataset %>%
  filter(test == "association" & mask == "unmasked" & session == "all" & control == "nature" & condition == "food") %>%
  unique() %>%
  group_by(process, subjectID) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
         processPE = paste0(process, "_ROI"),
         processPEV = paste0(process, "_PEV")) %>%
  ungroup() %>%
  select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map, process)) %>%
  spread(processPEV, dotProduct_std) %>%
  spread(processPE, meanProcessPEstd) %>%
  group_by(subjectID) %>%
  fill(contains("PEV"), .direction = "down") %>%
  fill(contains("PEV"), .direction = "up") %>%
  fill(contains("ROI"), .direction = "down") %>%
  fill(contains("ROI"), .direction = "up") %>%
  select(-c(roi, craving_ROI, craving_regulation_ROI, test, age, gender, mask)) %>%
  unique() %>%
  mutate(balance_ROI = cognitive_control_ROI - reward_ROI)

model_data_bmi = model_data %>%
  select(-contains("fat")) %>%
  na.omit()

model_data_fat = model_data %>%
  select(-contains("bmi")) %>%
  na.omit()

model_data_ema = model_data %>%
  select(-contains("fat"), -contains("bmi")) %>%
  right_join(., ema_tidy, by = "subjectID") %>%
  na.omit()

visualize

DVs

# source raincloud plot
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

# plot
model_data %>%
  left_join(., ema_tidy, by = "subjectID") %>%
  filter(!is.na(subjectID)) %>%
  select(subjectID, bmi, fat, enact_prop) %>%
  rename("BMI" = bmi,
         "body fat percentage" = fat,
         "enactment" = enact_prop) %>%
  gather(var, val, -subjectID) %>%
  ggplot(aes("", val, fill = var)) +
    geom_flat_violin(position = position_nudge(x = .1, y = 0), color = FALSE) +
    geom_boxplot(width = .1, outlier.shape = NA) +
    geom_point(position = position_jitter(width = .05), size = .5, alpha = .5) +
    facet_wrap(~var, scales = "free") +
    scale_fill_manual(values = pal_outcome) +
    scale_x_discrete(expand = c(.1, 0)) +
    coord_flip() +
    labs(y = "\nvalue", x = "") +
    theme_minimal(base_size = 14) +
    theme(legend.position = "none")

combined

combined_plot %>%
  filter(!process == "balance") %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature")),
         process = gsub("_", " ", process),
         test = factor(test, levels = c("neural signature", "ROI", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_wrap(~ test + process, nrow = 3) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

combined_plot %>%
  filter(!process == "balance") %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature")),
         process = gsub("_", " ", process),
         process = ifelse(!test == "ROI", sprintf("%s PEV", process), sprintf("%s ROI", process)),
         test = ifelse(!test %in% c("association", "uniformity"), "none", test),
         test = factor(test, levels = c("none", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = process)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~condition) +
  scale_fill_manual(name = "", values = pal_condition8) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

body composition plots

ROIs

Mean parameter estimates as a function of process, stimulus type, and control condition

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(condition, meanProcessPEstd, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(control, meanProcessPEstd, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

Correlations between variables and BMI or % fat

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(bmi, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(fat, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

PEVs

Mean correspondance as a function of process, stimulus type, control condition, and test

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(condition, dotProduct_std, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean correspondance\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(control, dotProduct_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean correspondance\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

Correlations between variables and BMI or % fat

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(bmi, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(fat, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

EMA plots

ROIs

model_data_ema %>%
  select(subjectID, contains("ROI"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
  unique() %>%
  gather(neuro_var, neuro_val, contains("ROI")) %>%
  mutate(neuro_var = factor(neuro_var, levels = c("balance_ROI", "cognitive_control_ROI",
                                                  "reward_ROI", "value_ROI"))) %>%
  gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
  group_by(ema_var) %>%
  do({
    plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
      geom_point(alpha = .5) +
      geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
      scale_color_manual(name = "", values = pal_condition) +
      labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
      theme_minimal(base_size = 14)
    
    print(plot)
    data.frame()
  })

PEVs

model_data_ema %>%
  select(subjectID, contains("PEV"), desire_pres0.prop, desire_str.m, goal_conf.m, resist.m, enact_prop, amount_Recoded.m) %>%
  unique() %>%
  gather(neuro_var, neuro_val, contains("PEV")) %>%
  mutate(neuro_var = factor(neuro_var, levels = c("cognitive_control_PEV", "reward_PEV", "value_PEV",
                                                  "craving_PEV", "craving_regulation_PEV"))) %>%
  gather(ema_var, ema_val, contains("m"), contains("prop")) %>%
  group_by(ema_var) %>%
  do({
    plot = ggplot(., aes(neuro_val, ema_val, color = neuro_var)) +
      geom_point(alpha = .5) +
      geom_smooth(method = "lm", alpha = .1, fullrange = TRUE) +
      scale_color_manual(name = "", values = pal_condition[2:6]) +
      labs(x = "\n mean parameter estimate", y = paste0(.$ema_var[[1]], "\n")) +
      theme_minimal(base_size = 14)
    
    print(plot)
    data.frame()
  })

correlations among ROIs and PEVs

dots_plot %>%
  filter(!grepl("craving", process)) %>%
  left_join(., betas_plot) %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(meanProcessPEstd, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nmean standardized parameter estimate", y = "mean standardized PEV\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

PEVs v ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(!grepl("dot", Var1) & !grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("roi", Var1) & grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

PEVs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("dot", Var1) & grepl("dot", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

run cross-validated models

descriptives

model_data %>%
  left_join(., model_data_ema) %>%
  gather(Variable, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop) %>%
  filter(!Variable == "balance_ROI") %>%
  group_by(Variable) %>%
  summarize(M = mean(val, na.rm = TRUE),
            SD = sd(val, na.rm = TRUE)) %>%
  filter(!is.na(Variable)) %>%
  mutate(Variable = gsub("_", " ", Variable),
         Variable = ifelse(Variable == "bmi", "BMI",
                    ifelse(Variable == "fat", "Body fat percentage",
                    ifelse(Variable == "enact prop", "Enactment", Hmisc::capitalize(Variable)))),
         Variable = factor(Variable, levels = c("BMI", "Body fat percentage", "Enactment",
                                                "Reward ROI", "Reward PEV",
                                                "Value ROI", "Value PEV",
                                                "Craving regulation PEV", "Craving PEV",
                                                "Cognitive control ROI", "Cognitive control PEV"))) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)

Variable	M	SD
BMI	23.19	3.36
Body fat percentage	25.43	8.00
Enactment	0.47	0.22
Reward ROI	0.06	0.50
Reward PEV	0.07	0.77
Value ROI	0.40	0.68
Value PEV	0.28	0.77
Craving regulation PEV	0.40	0.67
Craving PEV	-0.13	0.76
Cognitive control ROI	-0.13	0.50
Cognitive control PEV	-0.11	0.78

cor_mat = model_data %>%
  left_join(., model_data_ema) %>%
  select(bmi, fat, enact_prop, contains("ROI"), contains("PEV"), -balance_ROI) %>%
  gather(var, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop) %>%
  mutate(var = gsub("_", " ", var),
         var = ifelse(var == "bmi", "BMI",
               ifelse(var == "fat", "Body fat percentage",
               ifelse(var == "enact prop", "Enactment", Hmisc::capitalize(var)))),
         var = factor(var, levels = c("BMI", "Body fat percentage", "Enactment",
                                      "Reward ROI", "Reward PEV",
                                      "Value ROI", "Value PEV",
                                      "Craving regulation PEV", "Craving PEV",
                                      "Cognitive control ROI", "Cognitive control PEV"))) %>%
  arrange(var) %>%
  spread(var, val) %>%
  ungroup() %>%
  select(-subjectID) %>%
  as.matrix()

cors = Hmisc::rcorr(cor_mat)

cors$r %>%
  as.data.frame() %>%
  kable(format = "pandoc", digits = 2)

	BMI	Body fat percentage	Enactment	Reward ROI	Reward PEV	Value ROI	Value PEV	Craving regulation PEV	Craving PEV	Cognitive control ROI	Cognitive control PEV
BMI	1.00	0.66	-0.13	-0.12	-0.02	0.06	0.00	0.10	-0.02	-0.02	-0.06
Body fat percentage	0.66	1.00	-0.20	-0.07	-0.07	0.01	-0.04	0.02	-0.07	-0.03	-0.06
Enactment	-0.13	-0.20	1.00	0.21	0.31	0.23	0.29	0.03	0.21	0.01	0.05
Reward ROI	-0.12	-0.07	0.21	1.00	0.62	0.30	0.60	0.08	0.47	0.31	0.41
Reward PEV	-0.02	-0.07	0.31	0.62	1.00	0.58	0.91	0.19	0.81	0.47	0.59
Value ROI	0.06	0.01	0.23	0.30	0.58	1.00	0.65	0.28	0.51	0.06	0.16
Value PEV	0.00	-0.04	0.29	0.60	0.91	0.65	1.00	0.15	0.76	0.30	0.40
Craving regulation PEV	0.10	0.02	0.03	0.08	0.19	0.28	0.15	1.00	0.24	0.21	0.23
Craving PEV	-0.02	-0.07	0.21	0.47	0.81	0.51	0.76	0.24	1.00	0.45	0.54
Cognitive control ROI	-0.02	-0.03	0.01	0.31	0.47	0.06	0.30	0.21	0.45	1.00	0.91
Cognitive control PEV	-0.06	-0.06	0.05	0.41	0.59	0.16	0.40	0.23	0.54	0.91	1.00

cors$P %>%
  as.data.frame() %>%
  mutate_if(is.numeric, funs(ifelse(between(., .05, .1), "†",
                             ifelse(between(., .01, .05), "*", 
                             ifelse(between(., .001, .01), "**",
                             ifelse(. < .001, "***", "")))))) %>%
  kable(format = "pandoc", digits = 2)

BMI	Body fat percentage	Enactment	Reward ROI	Reward PEV	Value ROI	Value PEV	Craving regulation PEV	Craving PEV	Cognitive control ROI	Cognitive control PEV
NA	***		†
***	NA	†
	†	NA	*	**	*	**		*
†		*	NA	***	***	***		***	***	***
		**	***	NA	***	***	**	***	***	***
		*	***	***	NA	***	***	***		**
		**	***	***	***	NA	*	***	***	***
				**	***	*	NA	***	***	***
		*	***	***	***	***	***	NA	***	***
			***	***		***	***	***	NA	***
			***	***	**	***	***	***	***	NA

t-tests

model_data %>%
  select(-contains("z"), -balance_ROI) %>%
  left_join(., select(model_data_ema, subjectID, enact_prop)) %>%
  ungroup() %>%
  select_if(is.numeric) %>%
  map_df(~ broom::tidy(t.test(.x, mu = 0, na.action = "na.omit")), .id = "x") %>%
  rename("Variable" = x,
         "t" = statistic,
         "df" = parameter,
         "p" = p.value) %>%
  mutate(`Mdiff [95% CI]` = sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         Variable = gsub("enact_prop", "Enactment", Variable),
         Variable = gsub("fat", "Body fat percentage", Variable),
         Variable = gsub("bmi", "BMI", Variable),
         Variable = gsub("_", " ", Variable),
         Variable = factor(Variable, levels = c("BMI", "Body fat percentage", "Enactment",
                                                "reward ROI", "reward PEV",
                                                "value ROI", "value PEV",
                                                "craving regulation PEV", "craving PEV",
                                                "cognitive control ROI", "cognitive control PEV"))) %>%
  select(Variable, `Mdiff [95% CI]`, t, df, p) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)

Variable	Mdiff [95% CI]	t	df	p
BMI	23.19 [22.77, 23.60]	109.65	252	< .001
Body fat percentage	25.43 [24.44, 26.43]	50.28	249	< .001
Enactment	0.47 [0.42, 0.51]	20.58	96	< .001
reward ROI	0.06 [-0.00, 0.12]	1.90	261	.059
reward PEV	0.07 [-0.02, 0.17]	1.51	261	.132
value ROI	0.40 [0.31, 0.48]	9.38	261	< .001
value PEV	0.28 [0.18, 0.37]	5.84	261	< .001
craving regulation PEV	0.40 [0.31, 0.48]	9.52	261	< .001
craving PEV	-0.13 [-0.22, -0.04]	-2.74	261	.007
cognitive control ROI	-0.13 [-0.19, -0.07]	-4.14	261	< .001
cognitive control PEV	-0.11 [-0.20, -0.01]	-2.20	261	.029

specify model variables

options(na.action = "na.fail")
data_ctrl = caret::trainControl(method = "repeatedcv", number = 5, repeats = 5)

BMI

ROI models

bmi_roi_full = lm(bmi_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_bmi)
(bmi_roi_mods = MuMIn::dredge(bmi_roi_full))

MuMIn::get.models(bmi_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_bmi_roi = caret::train(bmi_z ~ reward_ROI + value_ROI,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

bmi_pev_full = lm(bmi_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_bmi)
(bmi_pev_mods = MuMIn::dredge(bmi_pev_full))

MuMIn::get.models(bmi_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_bmi_pev = caret::train(bmi_z ~ craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

bmi_null = lm(bmi_z ~ 1, data = model_data_bmi)
summary(bmi_null)

## 
## Call:
## lm(formula = bmi_z ~ 1, data = model_data_bmi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2100 -0.6542 -0.0558  0.4827  4.6115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01306    0.06326   0.206    0.837
## 
## Residual standard error: 1.006 on 252 degrees of freedom

bmi_combined = caret::train(bmi_z ~ reward_ROI + value_ROI + craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(bmi_aic = AIC(bmi_null, best_bmi_roi$finalModel, best_bmi_pev$finalModel, bmi_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*bmi_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

# ROI
summary(best_bmi_roi$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1355 -0.6606 -0.0899  0.4662  4.5052 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.03113    0.07236  -0.430   0.6674  
## reward_ROI  -0.30523    0.12990  -2.350   0.0196 *
## value_ROI    0.15583    0.09654   1.614   0.1077  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9976 on 250 degrees of freedom
## Multiple R-squared:  0.02502,    Adjusted R-squared:  0.01722 
## F-statistic: 3.207 on 2 and 250 DF,  p-value: 0.04214

best_bmi_roi

## Linear Regression 
## 
## 253 samples
##   2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 202, 201, 203, 203, 203, 203, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.9949716  0.03211436  0.7469028
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(best_bmi_roi$resample$Rsquared)

## [1] 0.03262441

sd(best_bmi_roi$resample$RMSE)

## [1] 0.1099513

# combined
summary(bmi_combined$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1994 -0.6654 -0.0886  0.5098  4.4311 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            -0.07322    0.07768  -0.942   0.3469  
## reward_ROI             -0.30708    0.12961  -2.369   0.0186 *
## value_ROI               0.11695    0.09989   1.171   0.2428  
## craving_regulation_PEV  0.14370    0.09792   1.467   0.1435  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9953 on 249 degrees of freedom
## Multiple R-squared:  0.03338,    Adjusted R-squared:  0.02173 
## F-statistic: 2.866 on 3 and 249 DF,  p-value: 0.03725

bmi_combined

## Linear Regression 
## 
## 253 samples
##   3 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 200, 203, 204, 203, 202, 204, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.9955154  0.03419623  0.7523059
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(bmi_combined$resample$Rsquared)

## [1] 0.03657507

sd(bmi_combined$resample$RMSE)

## [1] 0.1027421

% fat

ROI models

fat_roi_full = lm(fat_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_fat)
(fat_roi_mods = MuMIn::dredge(fat_roi_full))

MuMIn::get.models(fat_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_fat_roi = caret::train(fat_z ~ reward_ROI,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

fat_pev_full = lm(fat_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_fat)
(fat_pev_mods = MuMIn::dredge(fat_pev_full))

MuMIn::get.models(fat_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_fat_pev = caret::train(fat_z ~ craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

fat_null = lm(fat_z ~ 1, data = model_data_fat)

fat_combined = caret::train(fat_z ~ reward_ROI + craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(fat_aic = AIC(fat_null, best_fat_roi$finalModel, best_fat_pev$finalModel, fat_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*fat_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
    arrange(AIC))

summarize the overall best fitting

summary(fat_null)

## 
## Call:
## lm(formula = fat_z ~ 1, data = model_data_fat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9531 -0.5945  0.1396  0.6457  2.4201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01679    0.06321   0.266    0.791
## 
## Residual standard error: 0.9994 on 249 degrees of freedom

EMA enactment proportion

ROI models

enact_prop_roi_full = lm(enact_prop_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_ema)
(enact_prop_roi_mods = MuMIn::dredge(enact_prop_roi_full))

MuMIn::get.models(enact_prop_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_enact_prop_roi = caret::train(enact_prop_z ~ value_ROI,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

enact_prop_pev_full = lm(enact_prop_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_ema)
(enact_prop_pev_mods = MuMIn::dredge(enact_prop_pev_full))

MuMIn::get.models(enact_prop_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()

best_enact_prop_pev = caret::train(enact_prop_z ~ reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

enact_prop_null = lm(enact_prop_z ~ 1, data = model_data_ema)

enact_prop_combined = caret::train(enact_prop_z ~ value_ROI + reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(enact_aic = AIC(enact_prop_null, best_enact_prop_roi$finalModel, best_enact_prop_pev$finalModel, enact_prop_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*prop_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

summary(best_enact_prop_pev$finalModel)

## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.17754 -0.58154 -0.01067  0.53885  2.40418 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.09627    0.09985  -0.964  0.33745   
## reward_PEV   0.37158    0.11641   3.192  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9697 on 95 degrees of freedom
## Multiple R-squared:  0.09686,    Adjusted R-squared:  0.08735 
## F-statistic: 10.19 on 1 and 95 DF,  p-value: 0.001917

best_enact_prop_pev

## Linear Regression 
## 
## 97 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 77, 79, 77, 78, 77, 78, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.9722222  0.1498315  0.7677445
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE

sd(best_enact_prop_pev$resample$Rsquared)

## [1] 0.1591469

sd(best_enact_prop_pev$resample$RMSE)

## [1] 0.1075649

tables

bmi

bmi_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_bmi_roi$resample$Rsquared), 
                            mean(best_bmi_pev$resample$Rsquared), 
                            mean(bmi_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_bmi_roi$resample$Rsquared), 
                               sd(best_bmi_pev$resample$Rsquared), 
                               sd(bmi_combined$resample$Rsquared)),
                     RMSE = c(mean(best_bmi_roi$resample$RMSE), 
                              mean(best_bmi_pev$resample$RMSE), 
                              mean(bmi_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_bmi_roi$resample$RMSE), 
                                 sd(best_bmi_pev$resample$RMSE), 
                                 sd(bmi_combined$resample$RMSE)))

bmi_table = broom::tidy(best_bmi_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_bmi_roi$finalModel)$fstatistic[1], 
                             summary(best_bmi_roi$finalModel)$fstatistic[2], 
                             summary(best_bmi_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_bmi_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_bmi_pev$finalModel)$fstatistic[1], 
                             summary(best_bmi_pev$finalModel)$fstatistic[2], 
                             summary(best_bmi_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(bmi_combined$finalModel)$fstatistic[1], 
                             summary(bmi_combined$finalModel)$fstatistic[2], 
                             summary(bmi_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., bmi_aic) %>%
  left_join(., bmi_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "BMI") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

bmi_table %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Intercept	-0.07 [-0.23, 0.08]	0.08	-0.94	.347
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.37	.019
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Value ROI	0.12 [-0.08, 0.31]	0.10	1.17	.243
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Craving regulation PEV	0.14 [-0.05, 0.34]	0.10	1.47	.144
BMI	ROI*	721.74	0.03 (0.03)	0.99 (0.11)	Intercept	-0.03 [-0.17, 0.11]	0.07	-0.43	.667
BMI	ROI*	721.74	0.03 (0.03)	0.99 (0.11)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.35	.020
BMI	ROI*	721.74	0.03 (0.03)	0.99 (0.11)	Value ROI	0.16 [-0.03, 0.35]	0.10	1.61	.108
BMI	PEV	723.45	0.03 (0.03)	1.00 (0.13)	Intercept	-0.05 [-0.19, 0.10]	0.07	-0.67	.507
BMI	PEV	723.45	0.03 (0.03)	1.00 (0.13)	Craving regulation PEV	0.16 [-0.03, 0.34]	0.09	1.64	.102
BMI	Null	724.15	–	–	Intercept	0.01 [-0.11, 0.14]	0.06	0.21	.837

% fat

fat_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_fat_roi$resample$Rsquared), 
                            mean(best_fat_pev$resample$Rsquared), 
                            mean(fat_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_fat_roi$resample$Rsquared), 
                               sd(best_fat_pev$resample$Rsquared), 
                               sd(fat_combined$resample$Rsquared)),
                     RMSE = c(mean(best_fat_roi$resample$RMSE), 
                              mean(best_fat_pev$resample$RMSE), 
                              mean(fat_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_fat_roi$resample$RMSE), 
                                 sd(best_fat_pev$resample$RMSE), 
                                 sd(fat_combined$resample$RMSE)))

fat_table = broom::tidy(best_fat_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_fat_roi$finalModel)$fstatistic[1], 
                             summary(best_fat_roi$finalModel)$fstatistic[2], 
                             summary(best_fat_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_fat_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_fat_pev$finalModel)$fstatistic[1], 
                             summary(best_fat_pev$finalModel)$fstatistic[2], 
                             summary(best_fat_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(fat_combined$finalModel)$fstatistic[1], 
                             summary(fat_combined$finalModel)$fstatistic[2], 
                             summary(fat_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., fat_aic) %>%
  left_join(., fat_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Body fat") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

fat_table %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
Body fat	Null	712.18	–	–	Intercept	0.02 [-0.11, 0.14]	0.06	0.27	.791
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.07)	Intercept	0.02 [-0.10, 0.15]	0.06	0.38	.706
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.07)	Reward ROI	-0.14 [-0.39, 0.11]	0.13	-1.12	.265
Body fat	PEV	712.94	0.03 (0.03)	1.00 (0.08)	Intercept	0.00 [-0.12, 0.13]	0.06	0.05	.960
Body fat	PEV	712.94	0.03 (0.03)	1.00 (0.08)	Craving PEV	-0.09 [-0.26, 0.07]	0.08	-1.11	.268
Body fat	Combined	714.49	0.02 (0.03)	1.01 (0.06)	Intercept	0.01 [-0.12, 0.14]	0.07	0.19	.849
Body fat	Combined	714.49	0.02 (0.03)	1.01 (0.06)	Reward ROI	-0.09 [-0.38, 0.19]	0.14	-0.66	.508
Body fat	Combined	714.49	0.02 (0.03)	1.01 (0.06)	Craving PEV	-0.06 [-0.25, 0.13]	0.10	-0.65	.516

EMA enactment

enact_prop_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_enact_prop_roi$resample$Rsquared), 
                            mean(best_enact_prop_pev$resample$Rsquared), 
                            mean(enact_prop_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_enact_prop_roi$resample$Rsquared), 
                               sd(best_enact_prop_pev$resample$Rsquared), 
                               sd(enact_prop_combined$resample$Rsquared)),
                     RMSE = c(mean(best_enact_prop_roi$resample$RMSE), 
                              mean(best_enact_prop_pev$resample$RMSE), 
                              mean(enact_prop_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_enact_prop_roi$resample$RMSE), 
                                 sd(best_enact_prop_pev$resample$RMSE), 
                                 sd(enact_prop_combined$resample$RMSE)))

enact_table = broom::tidy(best_enact_prop_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_enact_prop_roi$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_enact_prop_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_enact_prop_pev$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(enact_prop_combined$finalModel)$fstatistic[1], 
                             summary(enact_prop_combined$finalModel)$fstatistic[2], 
                             summary(enact_prop_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., enact_aic) %>%
  left_join(., enact_prop_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Enactment") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

enact_table %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
Enactment	PEV*	273.28	0.15 (0.16)	0.97 (0.11)	Intercept	-0.10 [-0.29, 0.10]	0.10	-0.96	.337
Enactment	PEV*	273.28	0.15 (0.16)	0.97 (0.11)	Reward PEV	0.37 [0.14, 0.60]	0.12	3.19	.002
Enactment	Combined*	274.60	0.12 (0.13)	0.97 (0.13)	Intercept	-0.15 [-0.38, 0.09]	0.12	-1.25	.215
Enactment	Combined*	274.60	0.12 (0.13)	0.97 (0.13)	Value ROI	0.13 [-0.19, 0.46]	0.16	0.81	.419
Enactment	Combined*	274.60	0.12 (0.13)	0.97 (0.13)	Reward PEV	0.32 [0.05, 0.58]	0.14	2.33	.022
Enactment	ROI*	278.06	0.10 (0.09)	0.99 (0.14)	Intercept	-0.19 [-0.43, 0.05]	0.12	-1.59	.115
Enactment	ROI*	278.06	0.10 (0.09)	0.99 (0.14)	Value ROI	0.33 [0.04, 0.61]	0.14	2.27	.026
Enactment	Null	281.16	–	–	Intercept	-0.04 [-0.25, 0.16]	0.10	-0.42	.676

combined

bind_rows(bmi_table, fat_table, enact_table) %>%
  kable(format = "pandoc", digits = 2)

outcome	model	AIC	r2 (SD)	RMSE (SD)	term	b [95% CI]	SE	t	p
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Intercept	-0.07 [-0.23, 0.08]	0.08	-0.94	.347
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.37	.019
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Value ROI	0.12 [-0.08, 0.31]	0.10	1.17	.243
BMI	Combined*	721.56	0.03 (0.04)	1.00 (0.10)	Craving regulation PEV	0.14 [-0.05, 0.34]	0.10	1.47	.144
BMI	ROI*	721.74	0.03 (0.03)	0.99 (0.11)	Intercept	-0.03 [-0.17, 0.11]	0.07	-0.43	.667
BMI	ROI*	721.74	0.03 (0.03)	0.99 (0.11)	Reward ROI	-0.31 [-0.56, -0.05]	0.13	-2.35	.020
BMI	ROI*	721.74	0.03 (0.03)	0.99 (0.11)	Value ROI	0.16 [-0.03, 0.35]	0.10	1.61	.108
BMI	PEV	723.45	0.03 (0.03)	1.00 (0.13)	Intercept	-0.05 [-0.19, 0.10]	0.07	-0.67	.507
BMI	PEV	723.45	0.03 (0.03)	1.00 (0.13)	Craving regulation PEV	0.16 [-0.03, 0.34]	0.09	1.64	.102
BMI	Null	724.15	–	–	Intercept	0.01 [-0.11, 0.14]	0.06	0.21	.837
Body fat	Null	712.18	–	–	Intercept	0.02 [-0.11, 0.14]	0.06	0.27	.791
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.07)	Intercept	0.02 [-0.10, 0.15]	0.06	0.38	.706
Body fat	ROI	712.92	0.02 (0.02)	1.00 (0.07)	Reward ROI	-0.14 [-0.39, 0.11]	0.13	-1.12	.265
Body fat	PEV	712.94	0.03 (0.03)	1.00 (0.08)	Intercept	0.00 [-0.12, 0.13]	0.06	0.05	.960
Body fat	PEV	712.94	0.03 (0.03)	1.00 (0.08)	Craving PEV	-0.09 [-0.26, 0.07]	0.08	-1.11	.268
Body fat	Combined	714.49	0.02 (0.03)	1.01 (0.06)	Intercept	0.01 [-0.12, 0.14]	0.07	0.19	.849
Body fat	Combined	714.49	0.02 (0.03)	1.01 (0.06)	Reward ROI	-0.09 [-0.38, 0.19]	0.14	-0.66	.508
Body fat	Combined	714.49	0.02 (0.03)	1.01 (0.06)	Craving PEV	-0.06 [-0.25, 0.13]	0.10	-0.65	.516
Enactment	PEV*	273.28	0.15 (0.16)	0.97 (0.11)	Intercept	-0.10 [-0.29, 0.10]	0.10	-0.96	.337
Enactment	PEV*	273.28	0.15 (0.16)	0.97 (0.11)	Reward PEV	0.37 [0.14, 0.60]	0.12	3.19	.002
Enactment	Combined*	274.60	0.12 (0.13)	0.97 (0.13)	Intercept	-0.15 [-0.38, 0.09]	0.12	-1.25	.215
Enactment	Combined*	274.60	0.12 (0.13)	0.97 (0.13)	Value ROI	0.13 [-0.19, 0.46]	0.16	0.81	.419
Enactment	Combined*	274.60	0.12 (0.13)	0.97 (0.13)	Reward PEV	0.32 [0.05, 0.58]	0.14	2.33	.022
Enactment	ROI*	278.06	0.10 (0.09)	0.99 (0.14)	Intercept	-0.19 [-0.43, 0.05]	0.12	-1.59	.115
Enactment	ROI*	278.06	0.10 (0.09)	0.99 (0.14)	Value ROI	0.33 [0.04, 0.61]	0.14	2.27	.026
Enactment	Null	281.16	–	–	Intercept	-0.04 [-0.25, 0.16]	0.10	-0.42	.676

VIF

# BMI
car::vif(bmi_combined$finalModel)

##             reward_ROI              value_ROI craving_regulation_PEV 
##               1.101738               1.185011               1.085275

# %fat
car::vif(fat_combined$finalModel)

##  reward_ROI craving_PEV 
##    1.303633    1.303633

# enactment
car::vif(enact_prop_combined$finalModel)

##  value_ROI reward_PEV 
##   1.348297   1.348297

BMI	Body fat percentage	Enactment	Reward ROI	Reward PEV	Value ROI	Value PEV	Craving regulation PEV	Craving PEV	Cognitive control ROI	Cognitive control PEV
NA	***		†
***	NA	†
	†	NA	*	**	*	**		*
†		*	NA	***	***	***		***	***	***
		**	***	NA	***	***	**	***	***	***
		*	***	***	NA	***	***	***		**
		**	***	***	***	NA	*	***	***	***
				**	***	*	NA	***	***	***
		*	***	***	***	***	***	NA	***	***
			***	***		***	***	***	NA	***
			***	***	**	***	***	***	***	NA

BMI	Body fat percentage	Enactment	Reward ROI	Reward PEV	Value ROI	Value PEV	Craving regulation PEV	Craving PEV	Cognitive control ROI	Cognitive control PEV
NA	***		†
***	NA	†
	†	NA	*	**	*	**		*
†		*	NA	***	***	***		***	***	***
		**	***	NA	***	***	**	***	***	***
		*	***	***	NA	***	***	***		**
		**	***	***	***	NA	*	***	***	***
				**	***	*	NA	***	***	***
		*	***	***	***	***	***	NA	***	***
			***	***		***	***	***	NA	***
			***	***	**	***	***	***	***	NA

Visualization & Model Comparison

Dani Cosme & Rich Lopez

2020-03-06

load packages

define palettes

source data

standardize brain data

summarize outliers

merge data

tidy data for plotting

tidy data for modeling

visualize

DVs

combined

body composition plots

ROIs

PEVs

EMA plots

ROIs

PEVs

correlations among ROIs and PEVs

PEVs v ROIs

ROIs

PEVs

run cross-validated models

descriptives

t-tests

specify model variables

BMI

ROI models

PEV models

combined

summarize the overall best fitting

% fat

ROI models

PEV models

combined

summarize the overall best fitting

EMA enactment proportion

ROI models

PEV models

combined

summarize the overall best fitting

tables

bmi

% fat

EMA enactment

combined

VIF

BMI	Body fat percentage	Enactment	Reward ROI	Reward PEV	Value ROI	Value PEV	Craving regulation PEV	Craving PEV	Cognitive control ROI	Cognitive control PEV
NA	***		†
***	NA	†
	†	NA	*	**	*	**		*
†		*	NA	***	***	***		***	***	***
		**	***	NA	***	***	**	***	***	***
		*	***	***	NA	***	***	***		**
		**	***	***	***	NA	*	***	***	***
				**	***	*	NA	***	***	***
		*	***	***	***	***	***	NA	***	***
			***	***		***	***	***	NA	***
			***	***	**	***	***	***	***	NA